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SUMMARY 


The proper orthogonal decomposition technique (Lumley’s decomposition) is applied to the turbulent 
flow in a channel to extract coherent structures by decomposing the velocity field into characteristic eddies 
with random coefficients. 

In the homogeneous spatial directions, a generalization of the shot-noise expansion is used to deter- 
mine the characteristic eddies. In this expansion, the Fourier coefficients of the characteristic eddy cannot 
be obtained from second-order statistics. Three different techniques are used to determine the phases of 
these coefficients. They are based on (1) the bispectrum, (2) a spatial compactness requirement, and (3) a 
functional continuity argument. Results from these three techniques are found to be similar in most re- 
spects. The implications of these techniques and the shot-noise expansion are discussed in the appendix. 

The dominant eddy is found to contribute as much as 76% to the turbulent kinetic energy. In both 
two and three dimensions, the characteristic eddies consist of an ejection region straddled by streamwise 
vortices that leave the wall in the very short streamwise distance of about 100 wall units. 


1 INTRODUCTION 

The general recognition of the existence of organized motions or eddies in turbulent shear flows can 
be traced to the works of Theodorsen (1952) and Townsend (1956) over three decades ago. In the past 
20 years a great deal of insight into the characteristics of organized structures in turbulent shear flows 
has been gained, primarily by means of flow visualization and conditional- sampling techniques (Cantwell, 
1981). Some combined flow visualization and quantitative techniques (Kline, Reynolds, Schraub, and 
Runstadler, 1967; Falco, 1977) have demonstrated the significance of certain events (e.g., the bursting 
process) or structures in the turbulence production mechanism. 

Unfortunately, the present knowledge of organized motions has seldom been used in turbulence theo- 
ries or in quantitative models of turbulence. This is in part caused by the lack of a quantitative definition of 
organized structures and an objective means for assessing their contributions to turbulence stresses, partic- 
ularly their importance in the production of turbulence. In addition, most flow visualization studies have 
been carried out at low Reynolds numbers where the limited range of turbulence scales makes it easier to 
identify organized motions. Much of our knowledge of coherent motions is limited to those structures that 
can be seen in flow visualization experiments. It is desirable to have the means to extract coherent motions 
from fields and to evaluate their contribution to turbulence statistics, regardless of how chaotic the fields 
are. 


The need for quantitative descriptions of organized structures has led to the use of statistical tech- 
niques. One method, used by Townsend (1956), Grant (1958), and Perry and Chong (1982), and others, is 
to examine measured two-point correlation profiles for their consistency with a proposed structural model. 
Because of insufficient experimental data, the agreement of the models with all the components of the 
two-point correlation tensor has never been investigated. The data are usually deficient in the number of 
components of the tensor or the directions of probe separation. Extrapolation from this insufficient data 
can lead to confusion. For example, Moin and Kim (1985) have shown that the ability to infer structural 



information from the two-point correlation profiles is highly dependent on the direction of probe separa- 
tion. Casual inspection of two-point correlation profiles can also be misleading because one is seeking to 
extract information on the velocity vector from a tensor. 

The conditional or phase-averaging techniques (e.g., VITA technique of Blackwelder and Kaplan 
(1976)) are statistical methods designed to obtain the average structure that satisfies a prescribed condi- 
tion, usually at a single point. The difficulty with using conditional sampling methods is that the prescribed 
conditions are generally ad hoc and their relevance to actual flow conditions is unclear. Moreover, Adrian 
and Moin (1988) have shown that the two-point correlation tensor provides a good estimate of conditional 
velocity fields. Thus, in the neighborhood of the point at which the conditions are specified, the condition- 
ally averaged velocity can be extracted from the information contained in the correlation tensor. 

Most statistical techniques for extraction of organized structures from turbulent flows will produce a 
“structure” from virtually any stochastic field, whether or not structures of interest are in the field. Thus, 
the association of statistically derived structures with instantaneous events of dynamic importance must 
be predicated on independent knowledge that dynamically significant structures do exist. The result of 
any such statistical technique is an ensembl t-overciged structure or flow pattern. This flow pattern is of- 
ten confined to a small section of the flow domain with the surrounding structures averaged out, and the 
inherent symmetries in the statistics impose “artificial” symmetries on the resulting structures. Moreover, 
the interfaces of these structures are generally smeared when compared to the edges of instantaneous flow 
structures. Therefore, it is likely that averaged structures do not resemble the instantaneous flow structures 
in detail. The fundamental question is whether the structures that are deduced by statistical techniques are 
relevant This question will be addressed in section 6. Here we point out that for modeling purposes, one 
generally is not interested in every detail of the instantaneous structures so the ensemble-averaged structure 
may indeed be what is needed. 

In 1967 Lumley proposed a mathematically attractive definition of organized structures and a sta- 
tistical method for their extraction from stochastic turbulent velocity fields. The method is based on the 
decomposition of the fluctuating velocity field into a sum of mutually orthogonal eigenfunctions of the 
two point correlation tensor, weighted by random coefficients. The dominant (most energetic) eddy is de- 
fined to be the eigenfunction with the largest eigenvalue. The Karhunen-Loeve expansion (Loeve, 1955, 
Papoulis, 1965) is used for directions in which turbulence is statistically inhomogeneous. This decompo- 
sition has also been used for meteorological mapping (Obled and Creutin, 1986) and for data compaction 
and reduction (Ahmed and Rao, 1975). An important feature of this decomposition, that results from its or- 
thogonality properties, is that the contribution of the extracted eddies to second-order turbulence statistics 
can be determined. 

There has been skepticism that the long-time averaged, unconditioned two-point correlation tensor 
used in proper orthogonal decomposition can retain information about the highly intermittent unsteady 
structures which have been observed in turbulent wall layers. If a structure, no matter how intermittent, 
contributes a majority of the total integrated energy or Reynolds stress, then it will dominate the two-point 
correlation statistics; therefore, information about the structure will be retained in the correlation tensor. In 
fact, it has been shown (Adrian, Moin, and Moser, 1987) that linear estimates of classical conditional av- 
erages (quadrant II), which are computed from the two-point correlation tensor, are in excellent agreement 
with the actual conditional averages. Thus conditional averaging of this type yields little information that 
is not available from the two-point correlation tensor. However, if there are several dominant structures 
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with comparable energy, the situation is not clear, and the two-point correlation tensor may not provide 
adequate information about the structures. Of course, methods for extracting structures from the two-point 
correlation tensor also suffer from the same difficulties as do other statistical methods (e.g., smearing and 
artificial symmetries). Also, we will see in section 4 that characteristic eddy decomposition as formulated 
here is not unique and that external information must be supplied to uniquely determine the Fourier-phase 
coefficients of the resulting characteristic eddies. 

Although Lumley’s proposal was made 20 years ago, it has not been evaluated thoroughly because 
of a lack of the necessary experimental data; the complete two-point correlation tensor with at least one 
direction of probe separation is required. Payne (1966) and Bakewell and Lumley (1967) were the first to 
apply this technique. Payne (1966) used the two-point correlation measurements of Grant (1958) in the 
wake of a circular cylinder. Grant measured only the diagonal elements of the two point correlation tensor, 
Raa, (<» = 1,2,3) at three fixed positions. The remaining off-diagonal correlations were obtained by 
using the mixing-length assumption and the equation of continuity. The energy content of the dominant 
extracted eddy was not significantly larger than that of the next eigenfunction in the hierarchy. Because 
of anomalies in the results, particularly the presence of some negative eigenvalues (which represent the 
energy content of eddies), these results can not be considered conclusive. Bakewell and Lumley (1967) 
applied a simplified version of the decomposition theorem to obtain the most energetic eddy structure in 
the wall region ( y + < 40) of a turbulent pipe flow. The two-point correlation of the streamwise velocity 
component, Rn(r x ), was measured and decomposed. The other velocity components of the large eddy 
were obtained by using the mixing-length assumption and the equation of continuity. They reported that 
the largest eddy carries over 90% of the total streamwise turbulent intensity. Moin (1984), who performed 
the decomposition in only one and two dimensions, was the first to make use of the full correlation tensor. 
The correlation tensor was obtained from an underresolved numerical simulation of turbulent channel flow 
which made use of a turbulence model to account for the unresolved portion of the turbulence (large-eddy 
simulation). As in the previous studies, Moin found that the dominant eddy carried much of the turbulent 
kinetic energy (as much as 64%). Recently, Glauser, Leib, and George (1985) applied the scalar decom- 
position to the streamwise turbulent velocity fluctuations in an axisymmetric turbulent jet. Their results 
indicate that the dominant eigenfunction carries about 40% of the total streamwise turbulent intensity in- 
tegrated across the layer. 

In recent work by Herzog (1986), the correlation tensor R a p was measured in a pipe for a, /3 = 
1 and 3 ; the rest of the tensor was reconstructed from the continuity equation. These very ambitious 
experiments have produced the most comprehensive application of proper orthogonal decomposition to 
an experimental wall bounded flow, and they are the experiments most directly comparable to the current 
results. Herzog’s measurements were taken in a small subdomain near the wall; thus they are similar to 
the near-wall domain decomposition discussed in section 5. 

Possible causes for some quantitative differences between Herzog’s and the present results are outlined 
below. 

1 . Herzog measured the correlation tensor at only six points in the y and z directions and seven points 
in the x direction, this probably provides an inadequate number of degrees of freedom for the decompo- 
sition (see discussion at the end of section 3). Further, the spatial resolution of the measured correlation 
tensor varied greatly in the x and z directions, with good resolution for small separations (A x* = 19 and 
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A z* = 9) and progressively coarser resolution for larger separations (x* > 40 and z* > 20). This ne- 
cessitated the use of curve-fits to obtain representations of the correlation tensor which could be Fourier 
transformed. 

2. In Herzog’s experiment, the correlation was measured using multiple hot-film probes. To preserve 
the mathematical properties of the correlation tensor, these probes must have identical responses. Because 
they were not identical, some of the eigenvalues were negative; although this was not as serious as in the 
results of Payne (1966). The use of multiple probes can also lead to probe interference problems. Herzog 
minimized probe interference by using small probes with slender carriers, but it is still a concern, especially 
for small spatial separations. 

3. The correlations used for the current results were estimated from a smaller statistical sample than 
the sample used in the experiments, which could lead to some difference caused by statistical error. The 
impact of these difficulties is not clear. 

The current results and the results of Herzog are largely in agreement; for example, Aubry and 
Keefe (1987) found that the individual eigenfunctions they examined were similar in shape in most re- 
spects. However, there are some significant disagreements (see sections 3 and 5). The extent to which the 
disagreements are caused by the cited difficulties is a matter of speculation. Other potential causes of these 
disagreements are (1) the differences between a pipe and channel, and (2) the fact that Herzog’s pipe flow 
was not fully developed. 

The objective of the present work is to extract the characteristic eddies (as defined by Lumley ’s decom- 
position) in fully developed turbulent channel flow and measure their contribution to turbulence statistics. 
The characteristic eddies are those with maximal contribution to turbulent kinetic energy, but the theory 
does not maximize their contribution to the turbulence production mechanism (i.e., Reynolds shear stress). 
One of the results of this study is the contribution of the extracted eddies to the turbulent shear stress profile. 

A major difficulty with the decomposition technique is the treatment of the homogeneous directions, 
in which the Karhunen-Loeve decomposition is not useful. Lumley (1981) proposed using a generalization 
of the shot-noise decomposition (Rice, 1944). However, there is considerable arbitrariness in the speci- 
fication of this decomposition. In Lumley’s approach, the magnitudes of the Fourier coefficients of the 
decomposition are found easily, but the phases are more difficult. Lumley (1981) recommends using the 
third-order moments in the form of the bispectra to recover the phases. This proposal has not previously 
been implemented. One of the objectives of this study is to use additional statistical data to retrieve the 
phase information, and to compare the characteristics of the dominant structure to those obtained with other 
techniques. Some of the implications of using different specifications of the shot-noise decomposition have 
also been studied. 

The necessary statistical data is obtained from a data base generated by direct numerical simula- 
tion of turbulent channel flow (Kim, Moin, and Moser, 1987). This data base consists of instantaneous 
three-dimensional velocity and pressure fields collected at widely separated flow times. Calculations were 
performed at Reynolds number 3200 based on the centerline velocit y, U 0 , and channel half width, 8. The 
channel centerline corresponds to y* = yu T /u =180, where u T = J T w /p is the wall shear velocity. The 
computations were carried out with 128 x 129 x 128 grid points in the x, y, and z directions, respectively. 
The mean flow is in the x direction, and y is in the direction normal to the walls. Total averaging time was 
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about 1906 /t/o- The physical realism of the data has been verified by detailed comparison of statistical 
correlations and both instantaneous and conditionally averaged flow patterns with available experimental 
data. 


In section 2, the procedure for calculation of the two-point spectral density tensor is outlined. In 
section 3, the inhomogeneous turbulence decomposition or Karhunen-Loeve expansion in the direction 
normal to the walls, and its computational implementation, is presented. In section 4, the theoretical aspects 
of the shot-noise decomposition in the homogeneous directions are discussed. The characteristic eddy 
decomposition is applied in two and three dimensions in section 5, and the structure of resulting dominant 
eddies are examined, followed by conclusions and a general discussion, in section 6. 

We are grateful to Drs. Robert S. Rogallo and Sanjiva Lele for their comments on a draft of this 
manuscript 

2 CALCULATION OF TWO-POINT SPECTRAL DENSITY TENSOR 


Application of the orthogonal decomposition theorem to turbulent channel flow, with one direction 
of flow inhomogeneity and two homogeneous directions, requires knowledge of the two-point spectral 
density tensor. This tensor is calculated from the direct simulation data base described in section 1. The 
Karhunen-Loeve expansion requires the two-point velocity correlation tensor, 

Rij(r x ,y,y',r z ) =< Ui(x,y,z,t)uj(x + r x ,y',z + r z ,t) > (2-1) 

where u, ( t = 1,2,3) are the instantaneous turbulent velocity fluctuations in the streamwise (z), normal 
(y), and spanwise ( z ) directions, respectively. The < > denotes ensemble average which, because of flow 
homogeneity in x and z directions, is calculated by averaging in (x,z) planes as well as in time. It is 
actually more convenient to compute and use the two-point spectral density tensor k x ,y,y', k z ) , which 
is the Fourier transform of the two- point correlation tensor in r x and r t , that is 


<&ij(k x ,y,y',k z ) 


^2 If e~ ik%T *~ ik ’ T ‘ Rij{ r x ,y,y',r z ) dr x dr z 


where k x and k z are the wave numbers in the x and z directions. 


( 2 - 2 ) 


For computational purposes, the discrete Fourier transform of each instantaneous velocity field has 


been computed 

Ui(k Xl y,k z ,t n ) = J2ui(x,y,z,t n )e~' klX ~' k ‘* 
1,2 


( 2 - 3 ) 


The two-point spectral density is obtained from 


1 N ' 

Q>ij(k x ,y,y' ,k z ) = — u,( k x ,y, k z> t n )uj( k x ,y , k z ,t n ) 

N * «=i 

where N t is the number of instantaneous flow fields used for ensemble averaging and * denotes complex 
conjugate. Since d>, ; is the Fourier transform of , a real function, it is conjugate symmetric, 

%(kx,y,y',kt) ( 2 - 4 ) 
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The Navier-Stokes equations, and the boundary conditions in the channel flow are invariant with respect 
to two coordinate transformations, a reflection in the z direction (z mapped to — z and w becomes —w), 
and a reflection in the y direction about the centerline. For each velocity field in the ensemble, the velocity 
field obtained by any combination of these reflections is also included in the ensemble because they are 
equally valid solutions of the Navier-Stokes equations. This effectively quadruples the statistical sample 
and exactly enforces the following two symmetries in d>, ; (y = 0 at the centerline). 

%(k x ,y,y',k,) = ±<J> i; (fc x , -y, -y', k t ) (2-5a) 

®ij(k x ,y,y',k t ) = ±<&ij(k x ,y,y' ,-k z ) (2-5b) 

In equation 2-5a the minus sign is used when * = 2 or j = 2 but not when both are 2, while in equa- 
tion 2-5b, the minus sign is used when i = 3 or j = 3 but not when both are 3. The symmetry in 
equation 2-5b also results in a factor-of-two reduction in the computation required to perform the proper 
orthogonal decomposition in three dimensions Finally, <!>,, has the following symmetry because of its 
definition. 

®ij(k s ,y,y',k z ) = <& ; *(fc*,y',y, k z ) (2-6) 

These symmetries, and others which can be derived by combining them, will be used throughout the sec- 
tions that follow. 

3 THE KARHUNEN-LOEVE EXPANSION IN THE INHOMOGENEOUS 

DIRECTION 

A preliminary evaluation of the decomposition theorem can be performed by using the the decompo- 
sition of the correlation tensor, R i; (y,y'), for two points separated in the inhomogeneous direction only 
(y). This decomposition, known as the Karhunen-Loeve expansion (Loeve, 1955; Papoulis, 1965), is the 
foundation of the characteristic eddy decomposition for multiple dimensions presented in the next section. 
Note that while the decomposition may be performed in one or more dimensions, the underlying turbulent 
flow is always three dimensional and time dependent. The following material leading to equations 3-1 to 
3-7 can be found in Lumley (1970). It is presented here for continuity. 

Let v,( y) be a random vector function on a finite domain D. Given an ensemble of realizations of v,-, 
we wish to determine a deterministic vector function (or organized structure), </>,( y) , that has the highest 
possible mean square correlation with the members of the ensemble. That is, we wish to find </>,( y) that 
maximizes the ensemble average of the magnitude squared of the quantity 

_ fp v<(_ y) ( y) dy _ (3-i) 

a ~ (f D Mv)H(y)dy) l/2 ' 

Unless otherwise stated, in this paper the summation convention is implied for repeated indices. Note that 
in the above inner product only the shape, and not the magnitude, of (f>, is considered. It can be shown 
(Lumley, 1970) by the methods of calculus of variation that the desired <t> is a solution (eigenfunction) of 

J D Rij(y,y')<t>*}(y')dy' = Hdy) ( 3 * 2 ) 
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where Rij =< v i (y)vj(y') > is the two-point correlation function and <> denotes ensemble average. It 
can be shown that equation 3-2 does not have a unique solution; instead, there is a denumerable infinity of 
solutions, <f>\ n) (y ) , which can be normalized such that 


f ^ n) (y)^ ra) *(y)dy = 5nm < 3 ’ 3 ) 

J D 

Orthogonality implies that structures of different orders do not interact with each other in their contribution 
to second-order statistics. Each eigenfunction <f>\ n) ( y) is associated with a real positive eigenvalue X (n , 
and the eigenfunctions form a complete set. That is, the random vector field can be reconstructed from 


the eigenfunctions 

Vi(y) = £<v^ n) (y) 

n 


(3-4) 


where coefficients of different order are uncorrelated 


< °n°m ^ 



n= m 
nj m 


(3-5) 


Equation 3-4 is interpreted as the decomposition of the stochastic field v, into deterministic elements (or 
eddies) with random coefficients. It is expected that more deterministic velocity fields will be more effi- 
ciently represented by the expansion equation 3-4. An important consequence of equation 3-5 is that the 
contribution of each structure to the turbulent kinetic energy and turbulence stresses can be determined 


< v>( y) v,(y) >=£ A (n) 4>t* ( V) ^ ( V) < 3 ' 6) 

n 


(U1U . 

E = [ < ViVi > dy = £ ^ (n) ( 3 ‘7) 

where E is the total turbulence kinetic energy in the domain. The eigenvalues A (n) thus represent the 
contribution of each structure to total turbulent kinetic energy. 


It should be pointed out that the Karhunen-Loeve expansion can be formulated for any subdomain, 
Vi < y < y<i- In this case the limits of the integrals in equations 3-2, 3-3, and 3-7 are changed to j/i and 
y u and the eigenfunctions represent the characteristic structures in that subdomain. The division of the 
full domain of interest into two or more smaller regions may be advantageous for the convergence of the 
expansion and provide the means for further dissection of the flow field in a given region. For example, 
in turbulent boundary layers we may wish to find the characteristic structures in the wall and outer layers 
separately, rather than to search for one global structure for the entire flow. This is consistent with the 
general treatment of the problem of multiple scales in turbulent boundary layers. 

The above formalism is applied to the three-dimensional time dependent velocity field, u,(x, y,z t t), 
in turbulent channel flow. Our aim is to find the optimum representation, in the statistical sense outlined 
above, of the velocity field in the direction y normal to the walls. That is, given the velocity profiles u,( y) 
at all the ( x , z) locations and at all times, we seek deterministic functions that optimally represent the y 
variation of the velocity field. The desired s are obtained by substituting into equation 3-2 the two-point 
correlation tensor y') defined in equation 2-1 with r x = r, = 0. The integral equation 3-2 is solved 

numerically. The numerical approximation to the integral in equation 3-2 is given by 

r N 

/ /(y)dy^£ w ‘/‘ ( 3 * g ) 

Jo .=1 
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where /» is the value of / at a discrete grid point and u> is the weight function for the particular quadrature 
method used. In the present work we have used the trapezoidal rule with up to N = 129 nonuniformly 
spaced grid points for the entire domain, — 1 < y < 1 . The grid points are given by 


y, - -cos 


7r(; - 1) 


j - 1,2,... ,N 


(N - 1) 

The numerical approximation of the integral in equation 3-2 leads to an algebraic eigenvalue problem 

_»(») , __*(»») 

A (j) = X (n) <t> (3-9) 


where A is a 3 N x 3 N matrix and 

7 < "’ = [^(D.^OJ.^O) 

is the discretized nth eigenvector (of dimension 3 N), with 4>\ n) (i) the streamwise component of the the 
nth eigenfunction at the ith grid point. Because of the use of nonuniformly spaced grid points, A is not 
symmetric. However, a simple scaling transformation using the diagonal matrix with the diagonal, D = 
[ > • • ■ > \/un] transforms equation 3-9 into a symmetric eigenvalue problem. 

The numerical integration in equation 3-8 also could have been done using the Chebyshev polynomial 
representation that was used in the simulations of Kim, Moin, and Moser (1987), which would have been 
more accurate. The trapezoidal rule was selected because it allows much greater flexibility in solving the 
problem in subdomains. Use of the trapezoidal rule is justified because the correlation tensor is much 
smoother than the instantaneous velocity fields. 

It can be easily proven that the discrete system of equations in 3-9 preserves all the essential proper- 
ties of its continuous counterpart given by equation 3-2. All the eigenvalues are real and positive, and the 
eigenvectors are orthonormal when the discrete analog of equation 3-3 is formed by using the same quadra- 
ture rule that was used to approximate the integral in equation 3-2. In addition, the relations (eqs. 3.4-3.7) 
are also exactly satisfied when the discretized instantaneous velocity fields are represented in terms of the 
eigenvectors of A. These properties serve as a good check of the computer implementation of the method. 
It can also be shown for both the analytical and numerical problem that the eigenfunctions will satisfy the 
same boundary conditions as the velocity, and that in the three-dimensional decomposition (see section 4), 
the eigenfunctions will satisfy continuity (including dv/ dy = 0 at the wall). 

In this paper the eigenvalues will be arranged in descending order, with X (1) as the largest eigenvalue. 
We shall refer to <£ (1 > as the “dominant” eigenfunction or eddy. Whether <£ (1) is indeed dominant depends 
on the magnitude of X ( 1} relative to the magnitude of the remainder of the eigenvalue spectrum. This must 
be determined from the computations. For the one-dimensional formulation discussed in this section, the 
spanwise component of each eigenfunction, , is uncoupled from its streamwise and normal compo- 
nents. This is a reflection of the fact that for ( i = 1 , 2) , R u , R# are zero. Therefore, in the remainder of 
this section, denotes the spanwise component of the mth eigenfunction with zero streamwise and ver- 
tical components, and <f>\ l) and <$ denote the streamwise and vertical components of the Ith eigenfunction 
with zero spanwise component. However, when we refer to the total energy, E, it is the sum of all 3 N of 
the eigenvalues. In the multidimensional formulation discussed in the next section, all three components 
of each eigenfunction are fully coupled. 
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The first four eigenfunctions for the domain going from the wall to the centerline of the channel 
(0 < y* < 180) are shown in figure 1. Note that the eigenfunctions are normalized t o hav e a magnitude 
of unity, in accordance with equation 3-6; however, in figure 1, they are multiplied by y/SJ^ to allow com- 
parison of their relative contributions to turbulent stresses. The Karhunen-Loeve eigenfunctions generally 
behave in the same manner as other typical eigenfunctions; namely, the number of zero-crossings increase 
with the order of the eigenfunction. It is particularly significant that the streamwise, <f>\ ^ , and vertical, 
02 n) > components of the first three eigenfunctions have opposite signs throughout the domain and hence, 
make a positive contribution to turbulence production. This is not the case for some of the higher order 





Figure 1. First four one-dimensional eigenfunctions in the wall-to-centerline domain: , streamwise 

velocity (u); , normal velocity (v); , spanwise velocity (w). (a) First, (b) second, (c) third, and 

(d) fourth eigenfunction. 
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eigenfunctions. It is also interesting that <f>\ x \y) changes sign in the vicinity of the wall. A streamwise 
vortex located near the wall yields a spanwise velocity profile, w, which is similar to </4 1} . 

The contributions of the first three eigenfunctions to the total turbulent kinetic energy (eq. 3-7) in 
three separate domains are shown in table 1. The contribution of the dominant eigenfunction to turbulent 
shear stress, and hence to turbulence production, in each domain is indicated by the quantity P (1) /P (t) , 
where P (,) is the integral of the first term in equation 3-6 with t = 1 ,; = 2 , 

P<” = X (1) r^Wdy 

Jyi 

pO) is the integral of the total turbulence shear stress 

P (t) = f < ujU 2 > dy 

•'VI 


Table 1. Contributions of the one-dimensional eigenfunctions to energy and production. 


Domain 

n d 

A (1 >/£ 

\™/E 

\ {3) /E 

X (l)/ A ( 2) 

p(D/p(t) 

0 <y*< 40 

29 

0.61 

0.15 

0.08 

4.2 

1.03 

140 < y' < 180 

10 

0.44 

0.22 

0.20 

2.0 

1.94 

0 <y*< 180 

65 

0.32 

0.16 

0.08 

2.0 

0.66 


For the domain extending from the wall to the channel centerline, the dominant eigenfunction makes 
an appreciable contribution to the total turbulence kinetic energy, and its contribution to turbulence pro- 
duction is remarkably high. For the wall layer ( y' < 40), the contributions of the first eigenfunction to 
both turbulent kinetic energy and production are very significant. For all the cases tabulated, the dominant 
eigenfunction’s contribution to turbulence shear stress is significantly higher than its contribution to turbu- 
lent kinetic energy. Note that the first eigenfunction may contribute more than 100% of the Reynolds shear 
stress. While the formulation guarantees that convergence to the energy and to the turbulence intensities is 
monotonic, there is no such guarantee for the Reynolds shear stress. In fact, the decomposition emphasizes 
the shear stress of the lower eigenfunctions, requiring higher order eigenfunctions to contribute negatively 
to the Reynolds shear stress. In contrast, Herzog (1986) obtained Reynolds shear stress contributions which 
were less than 100% and did converge monotonically. The reason for this difference is not known (see the 
discussion of Herzog’s experiment in section 1). 

The convergence of the Karhunen-Loeve expansion for turbulence stresses as a function of the number 
of terms in equation 3-6 for the wall-to-centerline domain is shown in figure 2. Note that in this figure, and 
in the convergence plots shown in figure 3, the solid curve labeled “total” is the Reynolds stress taken from 
the direct numerical simulation of Kim, Moin, and Moser (1987). All velocities are nondimensionalized 
with the wall shear velocity, u T . The high Reynolds shear stress content of the lower order eigenfunctions 
is clearly evident. It appears that approximately 10 terms in the expansion are required to reproduce the 
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Figure 2. Convergence of - — , total eigenfunctions; , first eigenfunction; sum of first (a) 5, 

(b) 10, (c) 5, and (d) 5 eigenfunctions; sum of first (a) 10, (b) 20, (c) 10, and (d) 10 eigenfunctions. 

(a) u 2 , (b) v 2 , (c) w 2 , and (d) uv in the wall-to-centerline domain. 


turbulence stresses with reasonable accuracy. Note that for this case, the matrix A (eq. 3-9) is 195 x 195 
and possesses 195 orthonormal eigenvectors. 

The convergence of the expansion in the wall-layer ( y + < 40) is shown in figure 3. The convergence 
is remarkably fast; three to five terms in equation 3-6 are sufficient to reproduce all the turbulence stresses. 
In this case the matrix A is 87 x 87 . The convergence of any expansion in a subdomain is expected to 
be better than the convergence in the entire domain. However as shown in table 1, the K-L expansion 
converges faster in the wall layer than in a subdomain of the same size away from the wall, despite the 
fact that turbulence quantities vary most rapidly in the wall region. A possible explanation for better 
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Figure 3. Convergence of — , total eigenfunctions; , first eigenfunction; sum of first (a) 3, (b) 5, 

(c) 3, and (d) 3 eigenfunctions; sum of first (a) 5, (b) 10, (c) 5, and (d) 5 eigenfunctions, (a) u 2 , (b) v 2 , 

(c) w 2 , and (d) uv in the near-wall domain. 


convergence of the wall-layer eigenfunctions is that turbulence motions near a wall are more organized 
(deterministic) than outer layer turbulence resulting in a larger projection a in equation 3-1. 

The disparity of scales between the wall and outer layer suggests that better convergence may be 
obtained if the entire domain is split into two or more regions and the eigenfunctions for each region are 
calculated separately. To verify this assertion, the wall-to-centerline domain was split into the two parts, 
y+ < 40 and 40 < y* < 180 , and for each case, equation 3-2 was solved. The sum of the contributions of 
the dominant eigenfunction from each region is 45% of the total kinetic energy as compared to 32% when 
the domain was not split. No attempt was made to optimize the y location at which the domain was split. 
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Each velocity component has an alternative expansion that can be obtained by solving equation 3-2 
with only the autocorrelation of that velocity component. For example, streamwise eigenfunctions can be 
obtained from 

I Rn(y, j/) <p*( y') dy' = A y) (3-10) 

Jd 

Such an expansion for each velocity component clearly has a faster convergence rate than the one obtained 
from the full decomposition. For example, equation 3-10 was solved for the near-wall domain, and the 
dominant eigenfunction’s contribution is 74% of the total streamwise turbulent intensity. However, scalar 
decompositions for each component of the velocity fluctuations do not reveal any information about the 
contribution of the eigenfunctions to the Reynolds shear stress. 

Because the convergence of the Karhunen-Loeve expansion in representing energy is optimal, it is 
of interest to see how its convergence compares to that of the Chebyshev polynomials which are used in 
performing direct numerical simulations. To that end, the streamwise velocity component was decomposed 
in the full domain from one wall to the other by obtaining the eigenvalues of equation 3-10. These are 
compared with the energy carried by partial sums of the Chebyshev polynomials in the representation of 
the streamwise velocity fluctuations. The power of the Karhunen-Loeve expansion is evident in the first 
term; the first eigenfunction carries 23% of the energy; whereas the first Chebyshev polynomial carries 
only 4%. However, when one considers the number of terms required to represent the energy to a given 
tolerance, the performance difference is not as impressive. For example, to represent 90% of the energy, 
10 eigenfunctions and 12 Chebyshev polynomials are required, and to represent the energy to a part in 
10 3 , 35 eigenfunctions and 42 Chebyshev polynomials are needed. Thus the Karhunen-Loeve expansion 
is significantly advantageous if only one or two terms are to be retained; however, for accurate simulations 
of the type performed by Kim, Moin, and Moser (1987), the small improvement in accuracy obtained by 
using the Karhunen-Loeve expansion would not offset the increased computational cost such a scheme 
would entail. 

Finally, one should be cautious in drawing conclusions regarding the convergence of the expansions 
when too few grid points are used in equation 3-2 or when Rij is measured at only a few points y. This has 
often been the case in experimental measurements of two-point correlations (Herzog, 1986). Eigenvectors 
of equation 3-9 form a basis for a space of vectors of dimension 3 No, where No is the number of grid 
points in the domain. Thus, the number of terms in equation 3-6 required to recover all the turbulence 
kinetic energy at the grid points is always less than or equal to 3 No . One can be confident that a sufficient 
number of points have been used only if the number of terms required for convergence is significantly 
less than 3 No- Note that for scalar decompositions such as equation 3-10, the corresponding turbulent 
intensity at the grid points is recovered (by default) with less than or equal to No terms in the expansion. 
In table 1, the number of grid points in each domain is shown. 

4 THEORY OF CHARACTERISTIC EDDIES IN MULTIPLE DIMENSIONS 


The one-dimensional Karhunen-Loeve expansion described in section 3 provided some guidance to 
the significance of the “dominant” eigenfunction. The merit of this decomposition is evident in the wall 
layer where the dominant eigenfunction is indeed the major contributor to turbulence kinetic energy and 
production. However, eigenfunctions in one dimension do not represent eddies; and application of the 
decomposition method to the problem of identifying organized structures in turbulent flows requires that 
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the method be implemented in more than one dimension. In that way, the shape of the extracted eddy, as 
well as its contribution to turbulent stresses, can be determined. In this section we consider the theoretical 
foundations of the three-dimensional decomposition, with two spatial directions ( x and z) homogeneous. 
A two-dimensional decomposition can be developed similarly. 

For the three dimensional case, we wish to determine the eigenfunctions of the three-dimensional two- 
point correlation tensor R(r x ,y, y', r z ), where r s and r t are separations in the homogeneous directions 
1 and 2 . Since the Karhunen-Loeve eigenfunctions in homogeneous spatial directions are the Fourier 
functions (Lumley, 1981), we can equivalently consider the following eigenvalue problem: 

J D &ij(k s ,y,y\k t )4>*(k x ,y' t k,) dy' = \(k x ,k z )h(k x ,y,k z ) (4-1) 

where is the spectral density tensor discussed in section 2. The Karhunen-Loeve eigenfunctions are 
then £,-( k x , y, k z ) cxp(ik x x + ik z z) . These are not acceptable as characteristic eddies because the Fourier 
functions are not local in space, and we expect the eddies to be spatially compact. Moreover, the Fourier 
eigenfunctions are the eigenfunctions for dny statistically homogeneous system, so they do not reflect 
properties related to turbulence structure. The homogeneous spatial directions require a different treatment 
which, following Lumley (1981), will be based on a generalization of the shot-noise decomposition (Rice, 
1944). 

The eigenfunctions of equation 4- 1 have all the properties of the eigenfunctions developed in section 3 . 
In particular, eigenfunctions of different order are orthogonal and they can be normalized so that 

I dy = 6™ (4-2) 

and the Fourier transform of the velocity field can be reconstructed from the eigenfunctions with random 
uncorrelated coefficients 


and 


u%(k x ,y,k z ) — ^^Qnik x , k z )(f>\ ^ (.k x ,y,k z ) 

n 


< Qn( k x , k z ) o m ( k x , k z ) > — 


n {n) (k x> k z ), 

to, 


for n = m 
for n y m 


(4-3a) 


(4- 3 b) 


Note that the normalization condition equation 4-2 sets the magnitude of the complex eigenfunction <£, but 
leaves the phase unspecified; the phases may be set arbitrarily. 

To obtain the shot decomposition in the homogeneous directions, consider the first term in 
equation 4-3 for each k x and k z 

u\ l Hk x ,y,k t )=ai(k x ,k z )ti'\k x ,y,k z ) (4-4) 

which is the “dominant” term of the inhomogeneous decomposition, and has the maximum contribution to 
the kinetic energy at each wave number. Equation 4-4 can be inverse Fourier transformed to obtain 

u\'Hx,y,z) = j ^\x-x',y,z-z') a x{x',z , )dx'dz' (4-5) 
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The inverse transform of the coefficients oi is 01 which can be interpreted as a stochastic process in the 
homogeneous spatial directions, and the deterministic function ^ is the inverse transform of the first 
eigenfunctions for each wave number. This is the form of the generalized shot decomposition, in which 
the deterministic function ^ would be interpreted as a characteristic eddy which is distributed randomly 
in the homogeneous spatial directions (“sprinkled” or “scattered”) by the stochastic process ai . However, 
the Fourier transform of <j>\ X) obtained from equations 4-1 and 4-2 has arbitrary amplitude (set arbitrarily 
by the normalization condition in equation 4-2) and phase, and is therefore inappropriate as a characteristic 

eddy. 

To determine these quantities for the characteristic eddy (#), we seek a generalized shot decomposi- 
tion of u\ x \ 

u[ X) (x,y,z) = J (f> c i(x - x',y,z - z')g(x',z') dx'dz' (4-6) 

where the magnitudes and phases of the characteristic eddy are determined by some objective criteria, and 
g is the stochastic “sprinkling” process. The characteristic eddy ft can be related to </>, ( ) in equation 4-6, 
and g can be related to oi by the following relations: 

<f>i(x, y, z) = J (t>\ l) (x - x',y,z - z')f(x',z') dx'dz' (4-7a) 

and 

ai(x,z) = J f(x - x',z - z')g(x' ,z')dx'dz' (4-7b) 

where / is a deterministic function. Equation 4-7b can also be interpreted as a generalized shot decompo- 
sition of the homogeneous stochastic process oi . There are several subtleties to the decomposition equa- 
tion 4-6 and to the criteria by which the Fourier magnitudes and phases of # are determined. In particular, 
many criteria could be used to define the decomposition, so this decomposition is not unique. The criteria 
actually used to obtain the results presented in section 5 are discussed briefly in the following paragraphs; 
the appendix contains further discussion of this subject. 

To determine the magnitudes of the Fourier coefficients of / and therefore <ft, we require that g is 
“white,” in the sense that the integral of g in nonoverlapping intervals is uncorrelated (Lumley, 1981). 
This property determines the second-order statistics of the process g, which is assumed to have zero mean. 

< g(x,z)g(x',z') > = 8(x-x',z-z') (4-8a) 

and 

<g(k x ^,)g<ik' x ,k',)> = 6(k x -k' x ,k'-k' z ) (4' 8b) 

Obviously, g could be multiplied by any constant, which would only change the scaling of /. Other con- 
ditions on $ or / are possible (see the appendix); however, this is an appealing choice, since it makes the 
second-order statistics of g primitive. The function / then carries the second moment of ai , 

< a x (x,z)adx + 8x,z + 8z) >= J f(x,z)f(x + 8x,z + 8z) dxdz (4-9) 

Taking the Fourier transform and recalling equation 4-3b, the spectrum of / is obtained; 

\f(k x ,k,)\ 2 = \ iX) (k x ,k,) (4-10) 
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The Fourier magnitudes of / and therefore $ are thus determined; however, the phases of / cannot be 
determined from the second order statistics of g and 01 . 

Though the phases do not affect the contribution of the characteristic eddies to second-order statistics, 
they do determine the physical structure of the eddies. Results shown in section 5 use three different meth- 
ods to find the phases. The first method makes use of the third-order statistics of the stochastic processes 
oj and g, as suggested by Lumley (1981). In the second method, the characteristic eddy is required to 
be compact in space in a sense to be discussed below. The third method makes use of the fact that for 
continuous wave numbers (infinite computational domain), it is expected that the Fourier components of 
<f>' would be continuous functions of the wave numbers. These methods are discussed below. 

4.1 Phase determination from the bispectrum 

To make use of the third-order-statistics of oi and g, we follow Lumley (1981) and consider the 
bispectrum. The bispectrum of a stochastic process g, B g , is the Fourier transform of the three-point 
correlation function, R g (Lii, Rosenblatt, and Van Atta, 1976; Van Atta, 1979; Elgar and Guza, 1985). We 
have 


and 


Rg(r x ,r' x> r z ,r' t ) = < g(x, z)g(x + r x , z + r z )g(x + r' x , z + r' z ) > (4- 11a) 

B g (k x ,k' x ,k z ,k' z ) = <g{k I ,k z )g{k' x k' z )g\k x + k' x ,k z + k' z ) > ( 4 ' llb ) 


Using the Fourier transform of equation 4-7b and equation 4-10, the bispectrum of ai and that of g can be 
related 


B ai (k x , k' s , k z , k' z ) = B g (k x , k' x ,k z , k' z ) [A^( A;*, k z )X^ l \k x , k' z )\^\k x + k x , k z + k z )} 


1/2 


xe 


(4-12) 


where 6(k x ,k z ) is the phase of f(k x ,k z ) (i.e., f(k x ,k z ) = ^A (1) (fc*, We would like to 
require that B g be a real, positive constant. This is analogous to the “whiteness” property imposed on the 
second-order statistics of g, and implies that g has primitive third-order moments in the sense discussed 
above (see the appendix for further discussion). It is not possible to impose this condition, however, be- 
cause equation 4-12 completely determines the magnitude of B g for all wave numbers, and B g will not (in 
general) be constant. The condition that B g be real and positive requires that 

ip(k x ,k' x ,k z ,k' z ) = 6(k x ,k z ) + 6( k x , k z ) -6(k x + k' x , k z + k' z ) (4-13) 

where ip is the phase of B ai , which would allow the determination of the phases 6. However, in general, 
equation 4-13 has no solution since it represents order N 2 equations for the order N unknown 6' s, where N 
is the total number of Fourier modes. Instead we will require that equation 4-13 be satisfied approximately. 
This problem is encountered in optics, seismology, and signal analysis (Bartelt, Lohmann, and Wimitzer, 
1984; Matsuoka and Ulrych, 1984), and a variety of solution techniques have been proposed. For the results 
presented in section 5, equation 4-13 was solved in a weighted least-squares sense, where the weights are 
taken to be the magnitude of the bispectrum of a i (B ai )• However, for the two-dimensional decomposition, 
the phases were constrained to take on values of 0 or n (see below). 


16 



For two-dimensional decompositions (y and z), the symmetries of the flow result in a simplification 
of the bispectrum. The symmetry of dfy shown in equations 2-5 and 2-6 implies in the two dimensional 
case that 

%(y,y',k z )^±^(y,y',k z ) (4-14) 

where the plus and minus signs are as in equation 2-6. This implies that d\ ; is strictly real unless i = 3 
or j = 3 (but not both), in which case it is pure imaginary. This in turn implies that the eigenfunctions 
4>i(y, k z ) are strictly real for * = 1,2 and pure imaginary for i - 3. To compute the bispectrum of oi , the 
realizations of ai must be determined from u. This is done by using equation 4-3a and the orthogonality 
of the eigenfunctions to obtain 

oi (k z ) = J Ui(y,k z )fc(y,k t )dy (4-15) 

It is easily shown that when equation 4-15 is used to compute oi from the velocity field and oi from the 
z reflection of the velocity field (as discussed in section 2), the special properties of the eigenfunction just 
discussed imply that 

a,(k z ) =a\*(k z ) (4-16) 

The bispectrum of oi can be expressed 

B ai (k z ,k' z ) =< a\(k z )ai(k' z )a* x (k z + k' z ) > (4-17) 

The ensemble in this average includes velocity fields and their z reflection as discussed in section 2, which 
with equation 4-16 implies that the bispectrum B ai is real. Referring to equation 4-12, it is easily seen 
that if the phase of / is 0 or tt for all k z (i.e., / is real for all k z ), then the bispectrum of g will also 
be real for all k z and k' t , as desired. Also, equation 4-13 implies that if ^k z is added to any solution 
for 9(k z ), it will remain a solution. This simply represents a shift of the characteristic eddy in the z 
direction. Thus the problem of finding the phases 6{ k z ) of f is reduced to finding the sign of /. In the 
weighted least-squares solution of equation 4-13, the values of 6( k z ) are restricted to be 0 or tt, so that the 
bispectrum of g remains real. The least- squares problem is reduced to determining the sign of f(k z ), such 
that Y,k,k’ t B g { k z , k' z ) | B ai ( k z , k' z ) /B g ( k z , k! z ) | is maximum. This was solved by a discrete steepest-ascent 
algorithm. 


4.2 Phase determination from a compactness condition 

Given the magnitude of f as determined from equation 4-10, the phases can be determined by im- 
posing a condition on 4 rather than on g. In particular, we use the fact that when we speak of an eddy 
in a turbulent field, we mean a structure that is compact in space. This requirement of spatial compact- 
ness is one of the reasons the Fourier functions were rejected as characteristic eddies. The compactness 
of the characteristic eddy is very sensitive to the phases 0, so a compactness requirement can be used to 
determine the phases. This reasoning was employed by Herzog (1986). However, the compactness of a 
structure in space is difficult to define precisely. Instead, we consider some quantity of interest (e.g., one 
of the velocity components of the characteristic eddy) and determine the phases 6( k x , k z ) to maximize the 
maximum value of its integral over y. Such a requirement leads to a structure that is compact in space. To 
see why this is so, let I(x,z) be the integral of the quantity of interest in y. We define an integral area scale 
f I 2 (x,z) dxdz /J 2 (0,0), which is a measure of the area in the x and z directions over which I is signif- 
icantly nonzero. The integral in the numerator is independent of the phases of the Fourier coefficients of 
I, and is therefore independent of 0( k x , k z ) , if we consider quantities which are linear in the characteristic 
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eddy velocity. Therefore, when the magnitude of 1(0 ,0) is maximized, the area scale is minimized and 
the characteristic eddy is compact. 


It is easy to show that this compactness condition is satisfied when the Fourier coefficients of the 
integral in y of the quantity of interest are real and positive for all wave numbers. One could also define 
a compactness condition based on the value of a quantity of interest at a given y location. The integral 
condition stated above is used to avoid the need to specify the y location. In the results shown in section 5, 
the compactness conditions based on both the streamwise velocity u and the normal velocity v were used. 
In principle, the spanwise velocity w could also be used. However, the symmetries of 4>, ; (see sections 2 
and 4.1) imply that for k z = 0 , either <£ 3 or both 4>\ and <j >2 are identically zero. For all cases considered, 
^3 was zero, so a compactness condition based on w would have left the phases of the k z = 0 modes un- 
determined. For some higher eigenfunctions, it is 4>\ and <$2 which are zero, so a w compactness condition 
would be appropriate for them. 


4.3 Phase determination from wave-number continuity 


Finally, we expect that if the wave numbers are continuous, corresponding to an infinite domain in 
x and z, then ^>‘(k x , k z ) will be a continuous function of the wave numbers. Thus for the discrete wave 
numbers in the computation, we expect that ^(k x , k z ) will not be very different for neighboring wave 
numbers. This expectation can be enforced to determine the phases 6( k x , k z ). For ease of exposition, the 
procedure by which this is accomplished will be described for the two-dimensional (y and z) decompo- 
sition in which 6 depends only on k z . A similar though more complicated procedure is possible in three 
dimensions. The similarity of the eigenfunctions for neighboring wave numbers k z j and k z j + 1 results in 


the magnitude of the integral 

f (j>i(y,k Z j)4>?(y,k z j +i ) dy 


(4-18) 


being nearly unity. The phases 6( k z ) are set such that the integral equation 4-18 will be real and positive, 
which results in for neighboring wave numbers being as parallel as possible. This procedure leaves 
an overall phase undetermined. This phase is set (except for an overall sign) by the requirement that 
<£?( k t = 0) is real (i.e., 5(0) =0), which is necessary for to be real. The result of this specification of 
the phases will be called the zero-phase eddy because for all wave numbers, the phases of $ are the same 
in the sense just described. 


4.4 Overall algorithm 

A step-by-step description of the algorithm that was used to obtain the characteristic eddy given the 
two-point spectral density tensor is presented here. These steps are those required to implement the u 
compactness criteria; the other phase recovery techniques are similar. 

1. The matrix A (eq. 3-9) representing the numerical approximation of the integral in equation 4-1 
is formed and symmetrized for each wave number. 

2. The eigenvalues and eigenvectors of A are obtained. 
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3. For each wave number, the eigenvector associated with the largest eigenvalue is selected. 

4. The eigenvector is scaled so that its energy (computed as J D </>,<£* dy, with the integral approxi- 
mated by the trapezoidal rule) is equal to the associated eigenvalue. 

5. The eigenvector is multiplied by the complex quantity Y/lll where i = f D <t> i dy. Again, the 
integral is approximated by using the trapezoidal rule. 

6. When steps 1-5 are completed for each wave number, the eigenvectors are treated as the Fourier 
transform of a velocity field. The inverse Fourier transform in x and z is computed to obtain the physical 
space representation of the eddy. 

5 CHARACTERISTIC EDDIES IN MULTIPLE DIMENSIONS 

In this section, the decomposition discussed in sections 3 and 4 will be applied in both two and three 
dimensions. We consider the two-dimensional case separately because, as discussed in section 4.1, we 
will require the bispectrum, and only in the two-dimensional case do we have a statistical sample that is 
sufficient to obtain a reliable estimate of this quantity. In the two-dimensional case, we wish to extract the 
dominant eddy as viewed in planes perpendicular to both the wall and the flow direction ( y-z planes). The 
relevant laboratory experiment would be to illuminate a y-z plane with a light source in a smoke-filled 
channel flow and to identify characteristic eddies in this plane. For the two dimensional case, equation 4-1 
reduces to 

I &ij(y,y\k z )fa(y' ,k t ) dy = (5-1) 

Jy 

where the two-dimensional d>, ; is the Fourier transform of (eq. 2-1) with r x = 0. The eigenvalue 
problems (eqs. 4-1 and 5-1) are solved numerically for each wave number by using the method outlined 
in section 3. (Recall that the computation of the channel flow was done in a finite domain with a finite 
resolution, so there are a finite number of discrete wave numbers.) 


5.1 Contributions to second order moments 


The properties of the inhomogeneous decomposition developed in section 3 can be extended immedi- 
ately to the multiple-dimension case. In particular, the contributions of the eigenfunctions to second-order 
statistics can be determined (independent of the phases). 


< ti, U; > = £< u| n) u< n) >= EE (5-2a) 

n n k. 


and 


« = <**.*«> 

n k x k. 


(5 -2b) 


For the two-dimensional case, the sums with respect to k x and the dependence on k x are not present. 
In table 2, the contribution of the first three characteristic eddies to turbulent kinetic energy ( E (n) = 
k, k t ) ) are shown for the same domains as in section 3, for both two- and three-dimensional 

cases. Again, as an indicator of the contribution of the dominant eddy to turbulence production, the quantity 
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n ( 1 ) /P (t) is also tabulated, where 

Fl"> - EE dy (5-3) 

*, k, JVI 

and P (<) is defined in section 3. 


Table 2. Contributions of the two- and three-dimensional eigenfunctions to energy and production. 


Domain 

Dimensions 

E 0) /E 

e (2) /e 

e (3) /e 

e w /e {2) 

n (,) /p ( ‘) 

0 < y*< 40 

2 

0.64 

0.17 

0.07 

3.8 

1.05 

140 < y*< 180 

2 

0.46 

0.24 

0.18 

1.9 

1.55 

0 < y* < 180 

2 

0.38 

0.16 

0.09 

2.3 

0.77 

0 < y*< 40 

3 

0.76 

0.12 

0.06 

6.0 

1.04 

140 < y* < 180 

3 

0.57 

0.22 

0.12 

2.5 

1.53 

0 < y + < 180 

3 

0.50 

0.16 

0.08 

3.1 

0.75 


As in the one-dimensional case, the dominant eddy has a significant contribution to turbulent kinetic 
energy and production. The contribution of the dominant eddy to turbulent kinetic energy is somewhat 
higher than in the one-dimensional case, and is higher in three dimensions than in two dimensions. Note 
also that, as in the one-dimensional case, the energy contribution of the dominant eddy is much greater for 
the near-wall domain than for the same-size domain at the center of the channel. This again suggests that 
the near-wall turbulence is more organized than that away from the wall. 

The spectra of the first three eigenvalues for the wall layer (y< 40) and the domain extending from 
the wall to the centerline are shown in figure 4 for the two-dimensional case and in figures 5 and 6 for the 
three-dimensional case. In both two and three dimensions, the z spectrum of the dominant eigenvalue in 
the wall layer shows a pronounced peak at the nondimensional wave number k z & = 7.5 which corresponds 
to the wavelength \* t = 150 . The higher order eigenvalues have no such peak, and do not show a shift to 
higher k z values as would have been expected if small vertical scales were associated with small horizontal 
scales. The x spectra for the three-dimensional case fall off monotonically, similar to the streamwise 
turbulent kinetic energy spectra Kim, Moin, and Moser (1987). 

For the two-dimensional case, the contributions of the first two eigenfunctions to turbulent kinetic 
energy spectra, 

E (n) (k z> y) = f^(k z> y)f^(k z> y) = \^(k t )i^(y,k g )U^(y,k M ) (5-4) 

at selected y locations are shown in figures 7 and 8 for the near-wall and the wall-to-centerline domains, 
respectively. 

At most y locations the contribution of the dominant eddy is larger than the second eigenfunction at 
low wave numbers and is smaller at high wave numbers. An exception is shown in figure 8(a) where the the 


20 







Figure 6. Eigenvalue x spectra for the three-dimensional case. , first eigenvalue; , second eigen 

value; , third eigenvalue, (a) Near-wall domain and (b) wall-to-centerline domain. 


second characteristic eddy is in fact more energetic than the dominant eddy at the lowest wave numbers. 
The higher order eigenfunctions are associated with smaller scales in the y direction but not necessarily 
with smaller scale motions in the 2 direction. Note that the decomposition theorem requires that, at any 
wave number, the integral of E {V) over y must be larger than that of E (1) . These integrals are of course 
equal to respective eigenvalues plotted in figure 4. 

The spectra of the dominant eigenfunctions are broad banded, but as expected they are confined to 
a narrower range of energetic wave numbers than the velocity spectra (Kim, Moin, and Moser, 1987). 
Fourier decomposition of structures confined to a finite region in space leads to broad band spectra such 
as those shown in figures 7 and 8. Spectra with very pronounced peaks result from distinct Fourier modes 
which, in contrast to compact eddies, extend indefinitely in space. 

The spectra in figures 7(c) and 8(b) turn up slighdy at high wave numbers, while the spectra in fig- 
ure 8(c) are very noisy at high wave numbers. The noise in figure 8(c) is caused by roundoff error. The 
data used in these decompositions were stored with seven-digit accuracy, thus when spectra fall off by 
more than six or seven orders of magnitude, noise as in figure 8(c) is expected. The up-tums in figures 7(c) 
and 8(b) are caused by the finite resolution of the simulations. In the simulations, the dynamics of the 
highest wave numbers are not well represented because of the presence of the resolution cutoff. These 
anomalies in the spectra should have little or no effect on the simulations or the decompositions because 
they represent an insignificant amount of energy. 
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Figure 7. Contributions to the z energy spectrum of the first eigenfunction (□), and second eigenfunction 
(o) at (a) y* = 1 .5 , (b) y* = 10 , and (c) y‘ = 40 for the two-dimensional case in the near-wall domain 


(y* < 40). 






Figure 8. Contributions to the z energy spectrum of the first eigenfunction (□) and second eigenfunction 
(o) at (a) y * = 5 .4, (b) y* = 40, and (c) y* = 91 for the two-dimensional case in the wall-to-centerline 

domain (y* < 180). 



5.2 Characteristic eddies in the y-z plane 

We now turn to the problem of finding the two-dimensional characteristic eddy <f>^. The magnitudes 
of the Fourier coefficients of /, as defined in equation 4-7, are found from equation 4-10, which in two 

dimensions is reduced to 

| = (5-5) 

To obtain the phases of /, we consider the three techniques discussed in section 4: the bispectrum, com- 
pactness, and continuity in wave number space (zero-phase eddy). 

The bispectrum of ai (a real quantity) is plotted in figure 9(a) for the near- wall domain, the bispectra 



Figure 9. Contours of the bispectrum in the near-wall domain (y* < 40). k z and k' z are zero in the center 
of the plot, (a) the process oi (eq. 4-3). Negative contours are dashed. Contour increment is 10 4 . Note 
the twelve-fold symmetry of the bispectrum function. 
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Figure 9. Concluded, (b) the process g (eq. 4-6). Contour increment is 0.05. 


for the other domains are similar. In all three domains the bispectrum is positive over almost the entire 
wave number domain, there are only small negative regions along the axes where either k z , k z or k z + k t is 
zero. Because the bispectrum of oi is positive almost everywhere, the solution to the least- squares problem 
discussed in section 4.1 is not significantly different from the zero-phase eddy on which the sign of S Ql 
is based. In fact the only difference between the zero-phase eddy and the bispectrum-derived eddy is for 
the near-wall domain (0 < y* < 40 ), and it turns out that the difference occurs at only one wave number. 
With / defined by the approximate solution of equation 4-13, the bispectrum of g can be computed, and is 
shown in figure 9(b) for the near-wall domain. In all three cases, the bispectrum is again positive almost 
everywhere, as we require. 

Two compactness conditions have been used to determine the phases of / . They require that for 
j = l and 2 the integrals / #( y, z) dy obtain their maximum possible value somewhere in the domain; 
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in particular, we choose that they be maximized at z = 0 . These are the u and v compactness conditions, 
respectively. It is interesting that for the two domains including the wall ( y * < 40 and y* < 180), the u 
compactness condition is satisfied by the zero-phase eddy, although the v compactness condition is not. In 
the other domain (140 < y* < 180), the v compactness condition but not the u condition is satisfied by 
the zero-phase eddy. Note, that the u> compactness condition is not considered, as discussed in section 4.2. 

We now turn to the structure of the characteristic eddies in two dimensions. Figures 10 and 1 1 show the 
characteristic eddies with phases determined from the zero-phase condition. In the wall domain (y* < 40 ), 
the characteristic eddy consists of a narrow intense region of low streamwise velocity, which is about 50 
wall units wide. There is also a very narrow (25 wall-unit) region of intense velocity away from the wall. 
The velocity vectors projected into the y-z plane show that at the top of the domain there is a weak flow 
away from the center of the jet, making a vortex. 

In the wall-to-centerline domain (y* < 180) the characteristic eddy consists of low streamwise ve- 
locity with normal velocity away from the wall, similar to the eddy in the near- wall domain. Near the wall, 
the widths of these regions is the same as in the wall domain, and the velocity vectors reveal a pattern 




Figure 10. Two-dimensional zero-phase characteristic eddy in the near- wall domain ( y + < 40). (a) Con- 
tours of u and (b) velocity vectors projected into the y-z plane. Contour increment in (a) is 1. Negative 
contours are dashed. 
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similar to the near-wall domain. However, the center of the circulation region is somewhat higher than in 
the near-wall domain. Away from the wall, these regions become diffuse. 

In both the near- wall domain (y* < 40) and the wall-to-centerline domain, the characteristic eddy 
consists of an ejection and the accompanying low-speed region. If the sign of $ were reversed, the eddy 
would represent high-speed fluid moving toward the wall (sweep). The overall sign of the eddies is chosen 
for consistency with the bispectrum-derived eddy because the compactness conditions and the zero-phase 
condition leave the overall sign undetermined. The results of quadrant analysis in this flow indicate that 
ejections should be dominant far from the wall, while sweeps are dominant near the wall (Wallace, Eck- 
elmann, and Brodkey, 1972). Quadrant analysis indicates that the cross-over point for the dominance of 
sweeps and ejections is at y* py 12-15. The sign of the eddy is determined from third-order statistics, 
which also indicate the change in dominance from sweeps to ejections. In particular, the skewness of u 
goes from positive to negative at y* Pa 15 (Kim, Moin, and Moser, 1987). Thus, even in the near-wall 
domain (y* < 40 ), ejections are dominant for a majority of the domain, which explains the ejection char- 
acteristic eddy in the near-wall domain. To test these arguments, the characteristic eddy was computed in 
the domain y* < 10. The resulting eddy represented sweep motion as expected. 

Another interesting aspect of these characteristic eddies is that the magnitudes of the velocities are 
quite large. For example, in the near-wall domain the maximum (most negative) u velocity is — 16 (normal- 
ized by u T ), as compared to the maximum rms value of about 2.7. Note that the integral of the characteristic 
eddy velocity over z is a substantial fraction of the rms velocity fluctuation, and the characteristic eddy is 
confined to a relatively small region, resulting in the large local velocities observed. 



Figure 11. Two-dimensional zero-phase characteristic eddy in the wall-to-centerline domain (y + < 180). 
(a) Contours of u. Contour increment is 1. Negative contours are dashed. 
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Figure 11. Concluded, (b) velocity vectors projected into the y-z plane. 


The characteristic eddies just considered satisfied one of the two compactness conditions discussed 
before (the u condition) as well as the zero-phase condition. Characteristic eddies satisfying the v com- 
pactness condition and the bispectrum condition have also been generated and are nearly indistinguishable 
from figures 10 and 11. The main difference is that the maximum value attained by u is slightly smaller 
(9% in the near-wall domain) and the maximum value attained by v is slightly larger (6% in the near-wall 
domain). The reason for this similarity is that the differences between the criteria are either in the higher 
wave numbers ( k , > 50) where the spectrum has already fallen off significantly (see figure 4), or there is 
a difference at only one or two wave numbers. 

To evaluate the extent to which the characteristic eddy is observable, an instantaneous velocity field 
has been examined for structures. To find the locations of events corresponding to the characteristic eddy, 
sample functions of g( z) have been computed for the near-wall domain from the sample functions of oi 
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by using equation 4-7b. An example obtained at a particular x location and time is shown in figure 12(a). 
In this figure, there is a large, isolated positive excursion of g from zero at about z = 3.16, which should 
correspond to a structure similar to the characteristic eddy. A vector plot of the velocity in the y-z plane 
centered at this location is shown in figure 12(b) for comparison with figure 10(b). As in the characteristic 
eddy, a large streamwise vortex is to the left of the domain, with a flow away from the wall at the cen- 
ter of the domain. The shape of this vortex is somewhat different from the vortices in the characteristic 
eddy. The major difference between the instantaneous structure and the characteristic eddy is the lack of 
a counter-rotating partner in the instantaneous field. It is important to note that the appearance of pairs of 




Figure 12. (a) Sample functions of the g process (eq. 4-8) from the two-dimensional decomposition in the 
near-wall domain, (b) Velocity vectors projected into the y-z plane of the instantaneous velocity in the 
plane from which the sample function of g shown in (a) was taken. The plot domain is centered on the 
large peak in g appearing at 2 = 3.16. 


31 


counterrotating vortices in the characteristic eddy is a consequence of the symmetries in the statistics, and 
the techniques by which the phases were determined, and it does not imply that such pairs are prevalent 
in the instantaneous fields. Thus, the major difference in this case is caused by this “artificial” symmetry 
imposed on the characteristic eddy (see section 6 for further discussion of this point). 

5.3 Characteristic eddies in three dimensions 

In three dimensions, we must determine the function f(x,z) as defined in equation 4-7. When g 
satisfies equation 4-8, the magnitude of the Fourier coefficients of f are 

|/< *„*,)! = yfcHk^kz) (5-6) 

The phases of / were determined by the compactness conditions discussed in section 4 and the zero-phase 
condition. The bispectrum is not used in the three-dimensional case because the statistical sample available 
did not allow it to be accurately estimated. In the two-dimensional case, the compactness and bispectrum 
conditions produce the same results, but it is not clear whether this would occur in three dimensions. 

In addition to the zero-phase condition, the phases of / have been chosen to make the maximum 
values attained by f 4>\{x, y, z) dy as large as possible for i = 1 or * = 2 to obtain u- and v-compactness 
conditions. Unlike the two-dimensional case, the zero-phase characteristic eddy does not satisfy either 
of these compactness conditions. There is very little difference in the structure of the characteristic eddy 
between the two compactness criteria for the near- wall domain, but the magnitudes of the maximum values 
of the velocities do differ. Some discernible differences between the zero-phase eddy and the eddies derived 
from the compactness criteria in the wall-to-centerline domain will also be discussed. 

Contours of the streamwise and normal velocity components for the zero-phase eddy are shown in 
figures 13 and 14 for the near-wall ( y * < 40) domain. Shown are contours in the y-z plane at x = 0, the 
x-y plane at z — 0 , and the x-z plane at y* & 10 . The location x = 0 , z = 0 is where both the streamwise 
and normal velocities attain their maximum values, and can therefore be regarded as the center of the eddy. 
In the near-wall domain, the y-z contours reveal a region of low streamwise velocity approximately 50 
wall units wide that is coincident with a region of strong velocity away from the wall. This is consistent 
with the two-dimensional results. The x-y contours show that the low streamwise velocity region extends 
250 wall units upstream of x = 0 near the wall, decreasing to as little as 60 wall units upstream at y* = 40 . 
In contrast, the downstream extent of the low velocity region (150 wall units) does not decrease away from 
the wall. The region of intense vertical velocity has a streamwise extent of less than 100 wall units in the 
upstream and downstream directions. This is also evident in the x-z contours in the plane at y* = 10 . 

Contours for the u-compactness condition characteristic eddy in the wall-to-centerline domain are 
shown in figures 15 and 16. The v-compactness condition plots are similar, except for the magnitudes 
of the velocities. As in the two-dimensional case, the near-wall region of the characteristic eddy in the 
wall-to-centerline domain is remarkably similar to that in the near-wall domain. Contours in the y-z plane 
in figures 15 and 16 indicate that the region of low streamwise velocity and flow away from the wall are 
diffuse far from the wall, having a z extent as much as three times that near the wall. However, in the x-y 
plane, it is seen that these regions have smaller streamwise extents far from the wall than near the wall. Far 
from the wall, the low streamwise velocity region has an extent of about 100 wall units, as does the region 
of high vertical velocity. Contours of u and v in the x-z plane at y* = 95 reveal that upstream of x = 0 , 


32 




Figure 13. Contours of u for the three-dimensional zero-phase characteristic eddy in the near-wall domain 
(y* < 40). (a) y-z plane at x = 0 , (b) x-y plane at z = 0 , and (c) x-z plane at y* = 10 . Contour increment 
is (a) 7, (b) 7, (c) 5. Negative contours are dashed. 
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the highest values of u and v occur off the z = 0 axis. This is related to the pair of streamwise vortices in 
the characteristic eddy, because above these vortices, the maximum normal velocity does not occur on the 
midplane between them, but to the sides. The zero-phase characteristic eddy has also been computed; the 
major difference between it and the u-compactness eddy in figures 15 and 16 is the vertical location of the 
streamwise vortices (see below). 

Figure 17 contains plots of the velocity vectors projected into y-z planes located at x 4 = 0 ,± 18 and 
±36 for the zero-phase eddy in the near-wall domain. The plane at x = 0 reveals a pair of counter rotating 
vortices straddling the region of strong ejection from the wall. In the wall-to-centerline domain, a similar 
set of streamwise vortices occur (not shown). The centers of the vortices in the near- wall domain at x = 0 
are located at y 4 = 25 and are separated by 35 wall units. Upstream of x = 0 , the centers are closer to the 
wall, and downstream, they are farther away. At x 4 = 36 , the centers are beyond the edge of the domain. 

The y location of the vortex centers is plotted in figure 18. Note that the angle of inclination of the 
vortex increases with downstream location. At x* = —36, the angle is 10° while at x =18, the angle is 
31 °. Upstream of x 4 = -54 (plots not shown), the vortices are so weak that they are not discernible. Thus, 
in the wall region (y 4 < 40), the vortices have a streamwise extent of less than 100 wall units. This is 
consistent with the observations of streamwise vortices in instantaneous velocity fields made by Moser and 
Moin (1984), in which the vortices were observed to go beyond y 4 = 40 in a very short streamwise distance 
(150 wall units). In contrast, the streamwise vortices in the decomposition performed by Herzog (1986) 
were near the wall for a much longer streamwise distance (360 wall units), and had a nearly constant angle 



Figure 15. Contours of u for the three-dimensional u-compactness condition characteristic eddy in the 
wall-to-centerline domain (y 4 <180). (a) y-z plane at x = 0. Contour increment is 5. Negative contours 
are dashed. 
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of inclination of 5°. The reason for this discrepancy is not known, but the poor streamwise resolution of 
Herzog’s measured correlation tensor for moderate to large separations (see section 1) is a potential cause. 

Moser and Moin (1984) observed that there were many more solitary vortices than vortex pairs in 
instantaneous flow fields. It may appear that a contradiction exists between the counterrotating pair of 
vortices obtained by the decomposition and this prevalence of solitary vortices observed in instantaneous 
velocity fields. This is not the case. The technique outlined in sections 3 and 4 yields symmetric counterro- 
tating pairs of vortices, which are distributed through the homogeneous spatial directions by the sprinkling 
process g (see equation 4-6). As was observed in figure 12, there is no difficulty in generating a solitary 
vortex from a collection of vortex pairs when they are distributed by the appropriate g function. The fact 
that vortices are obtained by the current method indicates the importance (as measured by energy) of these 
vortices near the wall. Obtaining a pair of vortices, rather than a solitary vortex, calls for improvement in 
the phase recovery technique. 

Also shown in figure 18 are the vortex center positions for the characteristic eddies in the near- wall 
domain based on the u- and v-compactness conditions, as well as vortex center positions in the wall-to- 
centerline domain for all three phase determination methods. In the near-wall domain, the compactness 
conditions result in very similar curves, while the curve for the zero-phase eddy is somewhat lower. In 
the wall-to-centerline domain, the vortices are seen to be much farther from the wall than in the near-wall 
domain. At x = 0 , the vortex centers are at y * = 40 to 80, depending on the condition that was used to 
obtain the phase; while in the near-wall domain, the vortices were at y* = 20 to 30. The vortices resulting 



Figure 16. Contours of v for the three-dimensional u-compactness condition characteristic eddy in the 
wall-to-centerline domain ( y * < 180). (a) y-z plane at x = 0. Contour increment is 2. Negative contours 
are dashed. 
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Figure 16. Concluded, (b) x-y plane at z = 0, (c) x-z plane at y* = 10, and (d) x-z plane at y* = 95. 
Contour increment is (b-c) 2 and (d) 0.9. 
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from the compactness conditions are farther from the wall than the vortices for the zero-phase eddy. As 
in the near- wall domain, the vortices are inclined at about 10° upstream of x = 0. Downstream they 
arc inclined as steeply as 60°. In agreement with the near-wall results, the vortices are not discernible far 
upstream of x = 0 , and are near the wall for less than 100 wall units. Note that in both domains the vortices 
are far from the wall for x* > 30, but in figures 14(b) and 16(b), there is significant ejection near the wall 
as far downstream as x + = 180 . Thus, the far downstream ejection is below the streamwise vortex pair but 
is not straddled by them. Also, the inclined vortices observed here probably are related to, and certainly 
are consistent with, the hairpin eddies investigated by Moin and Kim (1985). 



Figure 17. Velocity vectors for the zero-phase eddy in the near-wall domain (y* < 40) projected into y-z 
planes at (a) x + = -36 , (b) x + = - 18 , and (c) x* = 0 . 
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Figure 18. Vertical (y) location of the center of the stream wise vortices as a function of streamwise dis- 
tance ( x ): •, zero-phase eddy; o, u-compactness eddy; □, v-compactness eddy. Upper curves are for the 
wall-to-centerline domain (y* < 180 ), lower curves are for the near-wall domain (y* < 40 ). 
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6 CONCLUSIONS 


The proper orthogonal decomposition technique has been applied to one-, two-, and three-dimensional 
decompositions of the turbulent channel flow. Velocity fields generated by direct numerical simulation 
(Kim, Moin, and Moser, 1987) were used to compute the full two-point velocity correlation tensor, 
Rij(r x , y, y',r z ), for this purpose. The technique was used to extract energetic organized structures from 
turbulence. Turbulent velocity fields were represented as a randomly weighted sum of the eigenfunctions 
of Rij. The resulting characteristic eddies were found to contribute as much as 76% to the kinetic energy, 
and even more to turbulence production. Thus the characteristic eddy indeed has special significance. 

Several different techniques were used to determine the Fourier phase coefficients of the characteristic 
eddy, including bispectral analysis and compactness conditions. The results were qualitatively similar for 
all techniques; however, quantitative differences, such as the y locations of the centers of streamwise 
vortices, were observed. 

In the two-dimensional case (spanwise and normal to the wall), the characteristic eddies consisted of 
a narrow ejection region accompanied by a region of low streamwise momentum. These were straddled by 
a pair of weak counterrotating vortices. It was found that characteristic eddies computed in the near-wall 
domain were very similar to the near-wall portions of the characteristic eddies that were computed in the 
full wall-to-centerline domain. 

Three-dimensional eddies were found to be similar to the two-dimensional eddies, in that they con- 
sisted of an ejection that was straddled by weak streamwise vortices. The low speed streak accompanying 
the ejection is about 400 wall units long and 50 wall units wide, while the region of strong vertical velocity 
is less than 200 wall units long. The counterrotating vortices are inclined at 10° near the wall and as much 
as 60° away from the wall, and they have a streamwise extent in the near- wall region (y* < 40) of less 
than 100 wall units. It was noted that while the streamwise vortices occurred in pairs, this did not imply 
that vortex pairs were in fact dominant in the instantaneous fields. The pairs resulted from symmetries in 
the statistics and from the techniques by which the phases were determined. In fact, as noted in section 5, 
solitary inclined vortices rather than pairs are most often observed near the wall in instantaneous velocity 
fields (Moser and Moin, 1984). 

The relationship of the characteristic eddies to structures observable in instantaneous flow fields is 
not clear, although it was demonstrated that the instantaneous structure associated with a large peak in 
the scattering function did have some qualitative similarity to the characteristic eddy. Presumably, the 
characteristic eddy represents some average of the most energetic events. It is clear however, that if the 
stochastic velocity field is composed of a single structure that is distributed in the homogeneous directions 
by a scattering process that is known, then the technique would extract that structure. Also, if the proper 
orthogonal decomposition eigenfunctions are indeed dominant, and if the scattering process is an impulse 
process (see the appendix) with impulses not too densely spaced, then the characteristic eddies would be 
visible in the flow. 

The structures in the instantaneous velocity fields are evolving in time. Unfortunately the three- 
dimensional characteristic eddy decomposition performed here does not recover these dynamics. An as- 
sessment of the role of the structures in the physics of near-wall turbulence requires some knowledge of 
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the evolution of the structures. The treatment of the dynamics of structures obtained by characteristic eddy 
decomposition is a topic of continuing research (Aubry, Holmes, Lumley, and Stone, 1988). 

However, whether the characteristic eddy is actually found in the instantaneous flow field, and if so at 
what frequency, may not be the important question. If the main objective of a study of coherent structures 
is to find a decomposition of turbulence into deterministic and stochastic portions, the characteristic eddy 
decomposition is certainly the most efficient method in the sense that energy and Reynolds stresses are 
reproduced with the fewest number of terms (“eddies")* In addition, governing equations for the coeffi- 
cients of the characteristic eddy eigenfunctions can be derived from the Navier-Stokes equations, which is 
a useful point for theoretical development. Other decompositions are also possible; for example, one could 
easily find the structures that are the most efficient contributors to the magnitude of vorticity fluctuations 
(Moser, 1988). 
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APPENDIX 

NOTES ON THE SHOT NOISE DECOMPOSITION 


In section 4 it was shown that several quantities had to be specified to fully determine $ in equa- 
tion 4-6. In particular, the second-order moments of the stochastic process g were specified (eq. 4-8). To 
determine the Fourier phases, we would have also liked to specify the third-order moments of g by re- 
quiring the bispectrum B g to be a real constant. In this appendix we examine the physical consequences 
of these specifications of the statistics of the process g. These statistical properties (equation 4-8 and B g , 
a real constant) are shared by many stochastic processes with vastly differing characters. In particular, 
processes which are independent on nonoverlapping intervals, satisfy these conditions. Two well-studied 
examples of such processes in one dimension (i) are the Gaussian white noise process and the Poisson im- 
pulse process. Gaussian white noise is a model for the velocity of a particle undergoing Brownian motion 
(note that B g is zero in this case). The Poisson impulse process is the g process used by Rice (1944) in the 
“shot-noise” decomposition. The Poisson impulse process can be written 

g(x) =^26(x - Xj) (A-l) 

(Papoulis, 1965), where x ; are random locations with uniform density along the x axis. More specifi- 
cally, the spacings between points are independent and exponentially distributed (this is necessary for the 
independent intervals property). The sample functions of these two processes have greatly differing char- 
acters. Thus the lower order statistical properties specified for g do not tell us very much about the sample 
functions of g. 

The Poisson impulse process is particularly interesting in connection with the current problem. If g 
in equation 4-8 is indeed an impulse process, then the stochastic process can be interpreted as the sum 
of characteristic eddies <f> < ' occurring at discrete random locations in the homogeneous directions. This 
is attractive because if the eddies are not too close together, they would actually be observable in the 
realizations of u,- !) , and perhaps in u, . It would thus be appropriate to require that g be an impulse process; 
that is, a process whose sample functions consist entirely of Dirac delta functions. An impulse process can 
be considerably more general than the Poisson process; in one dimension, a more general impulse process 

g( *) = 22 x ~ x i) ( A ' 2 ) 

i 

where the bj are not necessarily independent random variables, and the spacings between neighboring 
points x j are neither independent nor exponentially distributed. If it is arranged so that this process is 
homogeneous (as we require), its two-point correlation function would be given as 

< g(.x)g(x + r x ) >= v < b 2 > S(r x ) + S(r x )Rt,(r x ) (A-3) 

where v is the average density of points, S{ r x ) dr x is the probability that there is a point in the interval 
(r x ,r x + dr x ) given that there is point at 0; and R b ( r x ) is the correlation of the strength b of two impulses 
separated by a distance r x . As was discussed in the previous paragraph, the statistical conditions (equa- 
tion 4-8, and B g , a real constant) are certainly not sufficient to guarantee that g is an impulse process, and 
equation A-3 shows that they are also not necessary. Thus if our goal is to require that g be an impulse 
process, these specifications are not appropriate. 
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In general, it is not possible to require that g be an impulse process and satisfy equation 4-8, except 
in the limit of impulses which are dense in x. It may be possible however, to ask that g be an impulse 
process in some approximate sense. For example g in one dimension could be further decomposed to 
g(x) = gi(x) + g 2 (x) where g\ is an impulse process. This might be accomplished by requiring that 

^oi(x) - j f{x-x)9\(x)dxj ^ < A ~ 4 ) 

be a minimum. This was done by Chambers (1987) for a very restricted class of impulse processes g i . 
There are difficulties with this minimization problem because it should always be possible to reduce the 
error (eq. A-4) by increasing the average density of points in pi . However, allowing an arbitrarily large 
density of points is undesirable since the structures would then significantly overlap. It is also not clear 
how such a minimization might be accomplished in practice. 

Because we are concerned that g be as close as possible to an impulse process, it is interesting to see 
how well the sample functions of g that are determined from the conditions in section 4 approximate delta 
functions. To this end, some sample functions of g( z) from the two-dimensional case were computed from 
sample functions of oi by using equation 4-7b. An example is shown in figure 12(a). The sample function 
is certainly not made up of delta functions. However, it does exhibit regions of intense activity separated 
by regions of relative calm, and there are occasional, very large, short-duration excursions. This gives 
some indication that a decomposition of g into an impulse process and a remainder, as discussed before, 
might be a viable approach. 
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